function trans = markov2_noncausal_new_fast(image,T)
%second order markov feature extraction
%image=double(image);
dctresult1=image;
dctresult2=fliplr(flipud(image));
dctresult3=flipud(image);
dctresult4=fliplr(image);

array1x=dctresult1(:,1:end-1)-dctresult1(:,2:end);
array1y=dctresult1(1:end-1,:)-dctresult1(2:end,:);

array2x=dctresult2(:,1:end-1)-dctresult2(:,2:end);
array2y=dctresult2(1:end-1,:)-dctresult2(2:end,:);

array3x=dctresult3(:,1:end-1)-dctresult3(:,2:end);
array3y=dctresult3(1:end-1,:)-dctresult3(2:end,:);

array4x=dctresult4(:,1:end-1)-dctresult4(:,2:end);
array4y=dctresult4(1:end-1,:)-dctresult4(2:end,:);

array1x=round(array1x);
array1y=round(array1y);

array2x=round(array2x);
array2y=round(array2y);

array3x=round(array3x);
array3y=round(array3y);

array4x=round(array4x);
array4y=round(array4y);

array1x(array1x>T)=T;
array1x(array1x<-T)=-T;
array1y(array1y>T)=T;
array1y(array1y<-T)=-T;

array2x(array2x>T)=T;
array2x(array2x<-T)=-T;
array2y(array2y>T)=T;
array2y(array2y<-T)=-T;

array3x(array3x>T)=T;
array3x(array3x<-T)=-T;
array3y(array3y>T)=T;
array3y(array3y<-T)=-T;

array4x(array4x>T)=T;
array4x(array4x<-T)=-T;
array4y(array4y>T)=T;
array4y(array4y<-T)=-T;

trans = zeros(8, 2*T+1, 2*T+1, 2*T+1);
trans(1,:,:,:) = transmat_fast(array1x, T);
trans(2,:,:,:) = transmat_fast(array1y, T);
trans(3,:,:,:) = transmat_fast(array2x, T);
trans(4,:,:,:) = transmat_fast(array2y, T);
trans(5,:,:,:) = transmat_fast(array3x, T);
trans(6,:,:,:) = transmat_fast(array3y, T);
trans(7,:,:,:) = transmat_fast(array4x, T);
trans(8,:,:,:) = transmat_fast(array4y, T);
end

function v = transmat_fast(array, T)
  [M, N] = size(array);
  v = zeros(2*T+1, 2*T+1, 2*T+1);
  for i = 1:M-1
	for j = 1:N-1
	  v(array(i,j)+T+1, array(i,j+1)+T+1, array(i+1,j)+T+1) = v(array(i,j)+T+1, array(i,j+1)+T+1, array(i+1,j)+T+1)+1;
	end
  end
  s = sum(sum(v,3),2);
  for i = 1:2*T+1
    if (s(i) ~= 0)
	  v(i,:,:)=v(i,:,:) / s(i);
    end
  end
end